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Abstract 



Integral equation methods for the solution of partial differential equations, when coupled with suitable 
fast algorithms, yield geometrically flexible, asymptotically optimal and well-conditioned schemes in either 
interior or exterior domains. The practical application of these methods, however, requires the accurate 
evaluation of boundary integrals with singular, weakly singular or nearly singular kernels. Historically, these 
issues have been handled either by low-order product integration rules (computed semi-analytically) , by 
singularity subtraction/cancellation, by kernel regularization and asymptotic analysis, or by the construction 
of special purpose "generalized Gaussian quadrature" rules. In this paper, we present a systematic, high- 
order approach that works for any singularity (including hypersingular kernels) , based only on the assiimption 
that the field induced by the integral operator is locally smooth when restricted to either the interior or 
the exterior. Discontinuities in the field across the boundary are permitted. The scheme, denoted QBX 
(quadrature by expansion), is easy to implement and compatible with fast hierarchical algorithms such as 
the fast multipole method. We include accuracy tests for a variety of integral operators in two dimensions 
on smooth and corner domains. 
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1. Introduction 

One of the difficulties encountered in the practical application of integral equation methods lies in the 
need to evaluate integrals with singular or weakly singular kernels in complicated domains. For the sake of 
concreteness, we assume the computational task is to compute layer potentials such as the single and double 
layer potentials 



for target points x on a closed, smooth contour F C M^, whore G is the Green's function for an underlying 
elliptic PDE and fix' denotes the outward unit normal at x' . In the case of the double layer potential, it is 
typically the principal value of that is desired for x gT. 

In the present paper, we will restrict our attention to the Helmholtz equation 




(1) 



(2) 



A0-|-fcV = O, 



for which 



G{x,x') = ^H^^\k\x-x'\) 



(3) 
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where H^^'' denotes the Hankel function of the first kind of order 0. Hq^^ satisfies the Sommerfeld radiation 
condition 

lini (^-tk) = 0, 

where r = — and k € C with Im fc > 0. Using this Green's function, it is well-known that the kernels in 
([T]) , ^ are both logarithmically singular as operators acting on the boundary, and the quantities of interest 
are well-defined improper integrals. Off the boundary, the double layer potential must be treated with more 
care as the singularity is of the order l/|a; — a;'| and the limiting value has a jump of fJ.{x') at the point x' S F. 
In many applications, one would also like to compute the gradient of Sa or Dfi, which involves principal 
value and finite-part (hypersingular) integrals. 

When the target x is far from the boundary, the integrands in ([T]), ([2]) are smooth, and high-order 
quadratures can be obtained by standard methods. Difficulties are encountered only when x is either on 
or near the boundary. The problem of quadrature for singular or nearly singular integrals, of course, has a 



rich history and, we do not seek to review the literature here (see, for example, the texts Atkinson 1997 
Brebbia et al. 1984 Kress 1999| ). The most common approach is probably product integration, that is to 



say exact integration of the kernel multiplied by a piecewise polynomial approximation of the density a or 
fi on a piecewise smooth approximation of the boundary. Purely analytic rules, however, tend to be limited 
to a few singularities (such as log \x — x'\ in 2D or l/\x — x'\ in 3D) and low order approximations of the 



density and boundary. To handle the kernels HQ^\k\x — x'\) or e 



pik\x—x'\ 



analytic rules are often 



combined with numerical quadratures through the method of singularity subtraction. More precisely, it is 
easy to verify that 



K\k\x-x'\) 



2-K 



log \x 



and 



are smoother functions than the original kernels themselves and somewhat easier to integrate numerically. 
More generally, if the kernel Gi{x,x') can be integrated, say, on a flat surface by analytic means, then 
integrating 6*2 (x, x') — Gi{x, x') is an easier task if G2 and Gi have the same leading order singularity [Davis 



and Rabinowitzl 1984 Farina 20011 Johnson and Scott 1989 



Three other powerful approaches are (a) to design special purpose quadratures that integrate a specific 
class of singular functions with high-order accuracy 
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, (b) to find a change of variables that removes the principal singularity Bruno and Kunyansky 
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methods of type (b) are sometimes referred to as using singularity cancellation. In the complex analytic (or 
harmonic) case, we should also note that some remarkable interpolatory methods have been developed by 
Helsing and Ojala] |2008a for off-surface evaluation. 



The main purpose of the present paper is to introduce a rather different approach to the evaluation of layer 
potentials, based on the fact that the fields Sa or in ([T]), ([2| are locally smooth functions when restricted 
to either the interior or the exterior, although they may be discontinuous across the boundary. The scheme, 
denoted QBX (quadrature by expansion), is easy to implement, high order accurate, and requires only a 
smooth underlying quadrature scheme. This underlying smooth rule may be global or composite/panel-based 
and adaptive. The method is also compatible with fast hierarchical algorithms. QBX is an extension of the 



work of Barnett and Nguyen 2012| , who address the near but off-surface evaluation problem. 



The paper is organized as follows: In Section |2l we show how QBX follows naturally from considera- 
tions of potential theory. Section |3] describes the mathematical foundations of the method, and Section |4] 
demonstrates its numerical performance. A simple but complete description of the algorithm can be found 
in Section [331 Finally, we discuss additional details and potential extensions of the present work in Section 

13 
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(a) Error in the potential using the trapezoidal rule 
with 50 quadrature points. 



(b) Error in the potential using the trapezoidal rule 
with 100 quadrature points. 



Figure 1. 

functions. 



The potential Sa is computed using the trapezoidal rule, a simple, high-order quadrature for smooth 



2. Smooth, high-order quadrature 

Let us assume, for the moment, that we are given a smooth, simply connected closed curve F C M^, 
with a parametrization T = {'y{t) : < t < L}. We denote the interior of F by fi" and its exterior by 51+. 
We also assume that we have at our disposal an underlying quadrature rule capable of integrating smooth 
(non-singular) functions on F to high precision. In two dimensions, one option is the trapezoidal rule, since it 



is well-known to achieve superalgebraic convergence for smooth data on closed curves Davis and Rabinowitz 



1984 . 



A very natural question at this point is the following: for a target location x away from F, how well 
does the trapezoidal rule compute Sa{x) or D^{x)1 Certainly, the integrands in ([T]) and ^ are not actually 
singular in this situation, so the real question is how close x can be to F before accuracy is lost. Before 
analyzing this error more carefully, let us carry out a simple computational experiment for the curve 



lit) 



'I cos(t - 7r/4)(l + sin(2t)/2)\ 
sin(t-7r/4){l-hsin(2t)/2) ) ' 

with < t < 27r for a Helniholtz parameter k — 0.5. Using either 50 nodes (Fig. 1(a) I or 100 nodes (Fig. 



1(b) ), we plot the error in both the interior and exterior of F. The colors in these figures indicate the absolute 
value of the pointwise error for the Helmholtz single-layer potential Sa with a = \. 

These figures clearly suggest that the region in which the layer potential is inaccurate shrinks more or 
less in proportion to the grid spacing h (a fact well-known to practitioners of potential theory). To be a 
little more precise, let Ti^[Sa) denote the trapezoidal approximation of Sa using N points, and let 

E{x) = \Sa{x)-TN{Sa){x)\. 

For a fixed e, we define the "high-accuracy" region of the plane as the subset of where E{x) < e. This 
will, in essence, be all of with a tubular neighborhood of F removed. The extent of this tubular region 



depends on both h and e Barnett and Nguyen, 2012J. 



The fact that /i-refinement shrinks the region of inaccuracy, of course, is of no great value in evaluating 
layer potentials at points x on the curve F itself. For this, let us instead choose a point c off the surface with 

c = X + bhfij.1 , 

where hx is the unit normal to F at x. From our initial experiment, it is reasonable to expect that c is in 
the "high-accuracy" region. Assuming Sa is a smooth function in either the interior f2~ or the exterior r2+, 
it is easy to see that 



\TN{Sa){c) - Sa{x)\ = |(Tjv(5a)(c) - Sa{c)) + [Sa{c) - Sa{x))\ = 0(e + h) 



(4) 
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Figure 2. Geometric situation of Graf's addition theorem witli sources along the curve F, as used in (|5| and 
([6|. Note that x will reside on F further on in the discussion. 



since |c — x| = 0(h), under the assumption that the trapezoidal rule is accurate to precision e. In other 
words, the approximate value Tj^(S(j){c) is a first-order accurate approximation of the on-surface value 
Sa{x), within the error e. 

Remarkably, it is straightforward to improve matters even further. Instead of evaluating Sa{c), let us 
expand Sa about c to order p. The classical separation of variables representation of a smooth solution to 
the homogeneous Helmholtz equation in a disk centered at c takes the form 

oo 

cPix) = ^ aiMkp)e-''' (5) 

^ — — oo 

where (p, 9) denote the polar coordinates of the target x with respect to the expansion center c, and Ji is the 
Bessel function of order I (see Fig. [2|. For the single layer potential Sa, the coefhcients a; in the expansion 
([5]) can be computed analytically: 



Hl^\k\x' ~c\)e'^'^'a{x')dx', {I ^ -p, -p + I, . . . ,p) (6) 

r 



where {\x' — c\,6') denote the polar coordinates of the point x' with respect to c. These formulas follow 



immediately from Graf's addition theorem .Abramowitz and Stegun, 1965 



Hl^''\k\x^x'\)= J2 Hl^\k\x - c\)e'^''' Ji{k\x - c\)e-'^'^ , (7) 



— 00 



by interchanging the order of summation and integration. We note that Graf's addition theorem is generally 
applicable only if the target x is closer to the center than the source x': 

|a;-c| < |a;'-c|. (8) 

The integral defining ai is similar to that defining the original layer potential, except that HQ^\k\x' — c|) 

has been replaced with the more complicated but still smooth function Hi^\k\x' — c|)e''^ . Because of 
this smoothness, we evaluate a/ using the same trapezoidal rule T]^{ai). When seeking to evaluate Sa{x), 
however, we evaluate the local expansion ([5| instead. Figure 3(a) shows the result of using an expansion of 



order p = 3 superimposed on the naive use of the trapezoidal rule T/y to compute the single layer potential 
directly. The circular 'cut-out' regions in this and the following figures indicate where local expansions were 
used to approximate Sa. Note that the error in Sa computed by the local expansion is far smaller than 
the error of the naive computation throughout the circular region. In effect, the local expansion allows us 
to punch a disk-shaped hole into the tubular region of inaccuracy. This is precisely the idea underlying the 



close evaluation scheme of Barnett and Nguyen 2012 for targets near, but not on, F 



In this paper, we take the approach one step further. Namely, we investigate the use of the expansion 
([5| in evaluating the layer potential Sa{x) for x that actually lie on the curve F. 

Formally, it is worth noting that the radius of convergence of the local expansion about c is r = 
min^'gr 1^;' — c| as Graf's addition theorem requires it. Thus, we are seeking to evaluate a local expan- 
sion at its radius of convergence where the accuracy is most difficult to analyze. This difficulty, however. 
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(a.) p = 3, N : 



I quadrature nodes 



(c) p = 12, N = SO quadrature nodes 
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(b) p = 6, Af : 



I quadrature nodes 




(d) p = 12, A'' = 240 quadrature nodes 



Figure 3. The potential Sa is computed using the trapezoidal rule F with either A'' = 80 or A'' = 240 points, 
except in a disk of radius |c — a;| centered at an off-surface point c that lies in the "high-accuracy" region of the 
trapezoidal rule (here, approximately 3h away from the curve). Only a portion of the boundary F is plotted, and 
X is the point where the disk and F are tangent. We plot the error in the disk using various expansion orders p 
and numbers of quadrature nodes. 
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stems from the use of the addition theorem for a singular field (the potential due to a point source Hq^^). 
Figure 3(a) shows that the expansion is, in fact, accurate: it provides about four digits of precision uniformly. 
Loosely speaking, accuracy follows from the fact that the field induced by the layer potential is (one-sided) 
smooth in the interior and exterior domains or . The analytic issues here are somewhat involved and 
concern estimates on the decay of the coefficients ai in terms of the smoothness of the curve j(t) and the 
density a{t). Those estimates are established in Epstein et al. 2012 , and we will invoke them, as needed, 
below. 

If using an expansion of order p = 3 provides an accurate value for Sa{x), is it perhaps possible to 
obtain even more accuracy by further increasing p7 Figure |3(b)| shows the results of such an experiment. 
By setting p = 6, the accuracy of the potential in the vicinity of (and really also on F) increases from four 
to six digits. Further increasing the order to p = 12, however, causes a significant loss of accuracy (Fig. 
3(c) I. A consideration of the integrand in ^ shows why this occurs. As p is increased, both factors in the 



integrand Hi^\k\x' — c|)e''^ increase in complexity: e*'* by becoming more oscillatory, and Hi^\k\x' — c|) 
by becoming more sharply peaked. This combined effect leads to the resolution of the underlying trapezoidal 
rule being exceeded. Thus, in the experiment of Figure 3(c) the coefficients a/ are both large and wrong. 
Fortunately, once identified, this issue is easy to resolve. Indeed, simply increasing the number of points 
in the trapezoidal rule compensates for the added complexity of the integrands involved in computing the 
higher order coefficients, and with p = 12 more than ten digits of accuracy are achieved (Fig. |3(d)[ ). 

The sequence of experiments described thus far suggest a path to the high-order accurate evaluation 
of layer potentials as operators on the boundary. It also highlights one aspect of the scheme that requires 
careful analysis, namely the interplay between h and p. The local grid spacing h must be chosen small 
enough so that the coefficients in the local expansion ^ are computed with the necessary precision. We 
have concentrated in our experiments on a single boundary point x. To evaluate Sa everywhere on the 
boundary, we will simply introduce a large number of off-surface expansions centers whose corresponding 
disks cover a tubular neighborhood of F. We choose our expansion centers so that any desired target point 
is in the interior or on the boundary of one of these disks, enabling the application of QBX. If a target 
point is not in any of these disks, by definition it will be in the "high-accuracy" region associated with the 
original trapezoidal approximation. In the simplest approach, one may introduce an expansion center for 
each discretization node on the boundary. The procedure applied above to the single layer potential can be 
used just as well for the evaluation of integrals with hypersingular kernels. These are notoriously difficult 
using classical quadrature schemes. 

None of the observations made above change substantially if we replace the trapezoidal rule with another 
high-order quadrature. Figure |4] presents th e analog of Figure [T for a circle discretized using composite 



Gauss-Legendre quadrature. Using ideas from Barnett and Nguyen 2012 , we believe that the error contours 



(ignoring nodal oscillations) are the conformal images of the Bernstein ellipses Davis and Rabinowitz 1984 
associated with the integrand on each panel. 



QBX as a regularization scheme 

For readers familiar with multipole/partial wave expansions, the numerical results above may come as a 
surprise. After all, given a finite set of quadrature nodes (as in our computational examples), if the order of 
the local expansion were sufficiently high, it should converge to the field induced by a finite set of singular 
sources, namely the quadrature nodes. Instead, it is reproducing the continuous layer potential, even at the 
quadrature nodes themselves. There are two interpretations of this fact. 

For some, it is most natural to understand this in terms of series approximations of smooth functions, 
as introduced above. For others, it is perhaps useful to interchange the order of summation and integration 
and write 

Sa{x) ^ J Gp{x, x')a{x') dx' , 

where 

p 

Gpix,x') = Hl^\k\x' - c\)e'^^' Ji{k\x - c\)e-'^\ (9) 
i=-p 

for a target x on or near the boundary. That is, we can interpret the entire procedure as substituting the 
original Green's function G with Gp. It turns out that Gp is a surprisingly good filter. It regularizes the 
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(a) Error in potential from (smooth) composite 
Gauss-Legendre quadrature, with 5 panels consist- 
ing of 10 quadrature nodes each. 




(b) Error in potential from (smooth) composite 
Gauss-Legendre quadrature, with 10 panels consist- 
ing of 10 quadrature nodes each. 



Figure 4. The potential Sa computed using composite 10th order Gauss-Legendre quadrature. 



kernel in such a way that high-order accuracy is achieved without the need for additional correction. From 
this perspective, the need to decrease h with p is due to the fact that Gp itself is a more and more sharply 
peaked integrand as p increases. 



3. Mathematical Foundations of QBX 
3.1. Error analysis 

We turn now to the principal result justifying the use of QBX as a quadrature scheme. We restrict our 
attention to composite Gauss-Legendre quadrature, but the proof is analogous for any smooth high-order 
rule. In what follows, h will be used to denote the panel size in the composite Gauss-Legendre grid, rather 
than the point spacing used previously in discussing the trapezoidal rule. We apologize for this abuse of 
notation. 

Theorem 1. Suppose that T is a smooth, bounded curve embedded in M^, that Bj.{c) is the ball of radius r 
about c, and that Br{c)r\T = {x}. LetV be divided into M panels, each of length h and let q be a non-negative 
integer that defines the number ofnodes of the smooth Gaussian quadrature used to compute the coefficients 
a^^'^ according to the formula Id). For < /? < 1, there are constants Gp^p and C2q,p so that if a lies in 



the Holder space C^«''^(r), then 



p 



Sa{x) 



i=-p 



< 



4^) "I"" 



(10) 



Truncation error 



Quadrature error 



Proof. We begin by writing 



E 



Sa{x)- aiJi{k\x - c\)e''''- ) + E 



Ji{k\x - c\)e 



(11) 



The first term stems from using a truncated p-term expansion in Bessel functions to approximate Sa, while 
the second term is the error that stems from the numerical approximation of the coefficients in the truncated 
series. It is shown in Epstein et al. 20121 that the first error is of the order rP~^^ \\a\\cp.f>(r)- For the second 
term, we note that on a curve segnient 1'^ of length h, the standard estimate for g-point Gauss-Legendre 
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quadrature is Davis and Rabinowitz, 1984[ (2.7.12)] 



H'l^\k\x' -c\)e'^''' a{x')dx' 



< 



2q + 1 (2g) 

where D" denotes the n*'' derivative of the integrand with respect to the integration parameter along Ti . A 
straightforward combination of Stirling's approximation 

y2^n"+5e-" < n! < 20Fn"+5e-", 

summing over all panels, and bounds on the derivative allow us to write this term as 



p 

E 

i=-p 



aiJi{k\x ~ c|)( 



Ji{k\x - c\)e 



<C. 



9./3 



-) 

4r J 



2q 



Combining the two estimates yields the desired result. 

There are several aspects of the preceding theorem that are worth noting. 



(13) 



□ 



The two contributors to the error in the QBX approximation (10) are quite different. By placing an 



expansion center at a distance r = 0{h) away from the point x € T, the analytic truncation error 
is of order /i^"*"^. This error goes to zero under quadrature mesh refinement. In order for the second 
component of the error to be small, however, we need ^ < 1, so that r > h/A. A requirement of this 
type is essential, corresponding to the fact that if the expansion center is too close to the boundary 
relative to the discretization, accuracy will be lost. It is perhaps informative to set r = h/2, and write 
the error E from (11) in the form 

where e = (i)'^- Used in this manner, QBX is not classically convergent, but converges with controlled 
precision. 

• If one wants to achieve a classically convergent scheme, it suffices to refine r more slowly than h (say, 
with r = ^/h) . We prefer to keep the error components separate for the sake of clarity and because it 
permits additional tests of numerical consistency. 

Note that the estimate ( 10 ) explains the behavior of QBX discussed in Section [2j particularly the 



results shown in Fig. 3(c) 



• For the sake of simplicity. Theorem [T] assumes that the curve F is divided into equal-sized segments. 
In practice, with an adaptive discretization of the curve, a slightly different version of the result is 
needed. Since the difficulty in the error analysis is entirely local and the estimates are similar to those 
obtained above, we omit the rather cumbersome analysis. 

• In practice, one is often interested in evaluating the double-layer potential -D/x, or some derivative of 
Sa or Z?/i. Straightforward analysis shows that, for n derivatives of the Green's function, the error 
estimate in (10) is multiplied by a factor of r^". The use of QBX for such calculations is discussed in 
the next section. 

3.2. Derivatives, jumps, and principal value integrals 

Up to this point, we have focused on the calculation of the single-layer potential Sa. For the double layer 
-D/i defined in ([2]), the scheme is only slightly different. The coefficients in ^ are simply replaced by 



d 



H^^\k\x' - c\)e'^"' nix') dx'. {I = -p, + 1, . . . ,p) 



(14) 
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Because the scheme rehes on a local expansion of the potential, subsequent derivatives of Sa or Djj. with 
respect to the target location x are particularly easy to obtain by analytic differentiation of the local (Bessel) 
expansion. 

There is a complication which must be dealt with in evaluating operators other than the single layer 
potential (which is only weakly singular). QBX, by its construction, evaluates the one-sided limit of a layer 
potential, with the side determined by the location of the expansion center. In practice, however, one might 
want to compute the integral 



dG 

dfix' 



{x, x')fj,{x') dx' . 



(15) 



for X & As an operator acting on the boundary, D has a continuous kernel and D^i is well-defined Kress 



1999 . Dn is not, however, equal to its one-sided limit. Using the superscripts -I- and — to denote a point 



in the exterior or interior of F, respectively, the following jump relations are well-known Atkinson 1997 
Brebbia et al. 1984 Kress ~999| ). For the single-layer potential. 



Scf(x) 



lim Sa{x^), 



for its derivative, 



VxSa{x)\Y= lim Vx±Sa{x )±-na{x), 



where V xS<t[x) is defined in the principal value sense, and for the double-layer potential, 



(16) 
(17) 



D/i(a;)|r = lim D^{x )^-^{x). 

x^—^x Z 

For higher derivatives with respect to the target location a;, we have 

K 1 „ ^ 

dx,dxjS(j{x)\Y = lim d^±d^±Sa{x'^) ^ + '2hihj)a ± -{nii^ + ijhi) 



da 



(18) 



(19) 



where i is the unit tangent, assumed to satisfy the identity 



K, here is the curvature, 5 is the Kronecker symbol and s is arc length. The expression dxidxjS<T{x)\-c is 
defined in the Hadamard finite-part sense. 

Finally, in some settings, it is useful to consider derivatives of the double layer, tangentially oriented 
dipoles, and mixed source/target derivatives. For these, we have 



v{x') ■ Vx'G{x,x')cT{x')ds 



lim 



v{x') ■ \Ix'G{x^,x')<7{x') ds] T • v{x))a 



dx,D(T\T= lim {d^^±Da{x±))^\i~ 
x±^x i 2 ds 



(20) 



(21) 



In the former expression, v is the direction in which the source derivative is to be taken. When it is 
tangentially oriented, v{x') = t{x'), we denote the corresponding operator by R: 



R(t{x)^ / i{x') ■Wx'G{x,x')cr{x')ds. 



(22) 



The jump relations described in (19)-(21| are not so well-known (see, for example, [Kolm et al. 2003 



In summary, if the one-sided limit is the quantity of interest, then QBX computes that directly and no 
post-processing work is required. If, however, the principal value integrals Dfi{x), \7xSa{x), or the finite-part 
integrals in (21 1 are desired, then additional steps are required. The simplest scheme involves subtracting 
the relevant quantity from the QBX-derived one-sided limit. This retains the expected order of accuracy. 
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A second option is to compute both one-sided limits using QBX and average the quantities appropriately. 
That is, one can compute 

D^{x) = - [ lim D^{x+) + lim D^i{x~)] (23) 

by two applications of QBX. 

There are two drawbacks to the latter approach and one advantage. First, it makes the scheme approx- 
imately twice as expensive as using a one-sided limit. Second, the limiting values obtained from the two 
sides of r can vary noticeably in their accuracy for a given choice of smooth rule and expansion order. As a 
result, the error in the principal value computed by averaging is dominated by the worse of the two limits. 
The advantage of using the two-sided limit is that the Nystrom approximation of the operator Dfi{x) is 
much better behaved spectrally. We discuss this issue in some detail in Section [33j When solving integral 
equations, we believe this advantage outweighs the other considerations. 

3. 3. Informal description of the algorithm 

This section provides a complete description of the steps required to implement QBX. We assume that 
we are given a smooth curve T subdivided into M panels Fi, . . . , F^f of arc lengths hi, ... , Hm, respectively 
and that n denotes the outward normal to F. 



Set up parameters 

1. Fix the desired accuracy e. 

2. Choose local expansion order p (so that Sa will be computed to the order of accuracy p+1). 



3. Choose q and r such that ( 10 ) is approximately satisfied to precision e (assuming the un- 



derlying smooth rule is composite Gauss- Legendre quadrature). For points on panel rn, a 



value of Tm = hm/'2 works well in practice. (See Section 3.1 for details.) 

Compute one-sided limit 

4. For each target point xj G F„j C F: 

(a) Fix the expansion center Cj := Xj ^ {h„i/2)n, with (— ) corresponding to seeking the 
interior limit and (-f ) corresponding to seeking the exterior limit. 

{ If Cj is too close to any other panel n ^ m, refine the quadrature ("source") grid 
(thereby shrinking h,n o,nd moving Cj closer to Fj until this is no longer case. } 

(b) Compute the expansion coefhcients. For example, for the single layer potential, 

a,,i J^Hl'\k\x'-c,\)e^'''a{x')dx' 

for I = —p, . . . ,p using the underlying qth order accurate rule. (See Fig. [2] for the 
definitions of 9,0' .) 

(c) Evaluate the local expansion at c^: 

p 

Uj := ^ aj^iJi{k\x - c\)e~'''^. 
i=-p 

5. If the desired integral is a principal value or finite-part integral that has a jump condition. 



use the appropriate expression from Section 3.2 to subtract the appropriate term from the 



one-sided limit (or repeat the calculation with a center on the opposite side and average the 



two sided limits as in (23l) 



A few observations are in order: 



• In practice, the error from QBX is greater when the expansion center for a target point x lies on the 
concave side of the curve rather than the convex side. (See Fig 



and Nguyen 2012 for analytic insight.) 



1(a) for an illustration and Barnett 
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The algorithm contains a few nested loops, allowing for algorithmic variation. One can save storage, 
for example, by avoiding the allocation of memory to the expansion coefficients aj^i. Each "source" 
point can compute its contribution to Uj directly. 

The above algorithm directly applies the layer potential operator to a given density a{x'). It is 
straightforward to modify the algorithm to compute and store all (or near neighbor) quadrature weights 
in a table, hence constructing the Nystrom matrix approximating the integral operator (see Section 



3.4.2) 



3.4- Remarks on grids 

The QBX procedure does not require tight coupling between the "source" and "target" grids. By "source" 
grid, we mean the set of points on F where the density and Green's function are sampled in computing the 
local expansion coefficients using the underlying smooth quadrature rule. By "target" grid, we mean the set 
of points along F where we seek the value of the layer potential. 

3.4.1. Adaptive boundary grids 

In our discussion thus far, we have implicitly assumed the subintervals used to divide the boundary F 
are all of approximately the same length. Many applications, of course, are best addressed using some form 
of adaptive mesh refinement along the curve to resolve either complicated data or to discretize an unknown 
but complicated density. 

The main issue for the application of QBX with such grids concerns the location of the expansion centers 
and the validity of the error estimate in ( [To| . Fig. [S] illustrates the issue, under the assumption that centers 
are chosen using the rule r « /ii/2, /12/2, /13/2 on three successive panels Ii, 12,13. Consider, now, the 
expansion center indicated by c, which is at a distance r = ft.2/2 from F. The filled triangles in the figure 
illustrate the angles spanned by adjacent source quadrature nodes. The resolution provided in evaluating 
expansion coefficients at c by the sources on the panel Ii is clearly much lower than that provided by the 



sources on I2 or J3. Moreover, the assumption that r > /ii/4, which is essential in (10) in order for the error 
to be small, is violated. 

Fortunately, this problem is straightforward to address: one simply requires sufhcient sampling on Ii for 
the error estimate to guarantee high precision. There are several possible strategies in terms of implementa- 
tion, and we list two here. 

• If adjacent panels have substantially different lengths, interpolate the source density on the larger one 
to a fine grid that matches the resolution of the smaller one on the fly. 

• In discretizing the boundary, require that no two adjacent panels differ in length by more than a factor 
of two and increase the number of points q per panel by a factor of 2. 

We use the second (simpler) strategy for the data presented in Section [3] 
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3-4.2. Grids for solving integral equations 

When solving integral equations with a Nystrom method, a common grid is used for both sampling the 
unknown density and evaluating the resulting layer potential. This coincides with what we have referred to 
as the target grid, which should resolve the curve and the density to the desired precision. This grid may not 



be sufficiently fine to satisfy the requirements (10) and ([8|). Under those conditions, in the QBX procedure, 
one simply needs to interpolate the density from to a finer grid, which becomes what we have referred to as 
the source grid. 

3. 5. Spectral structure of operators approximated by QBX 

In this section, we consider the spectral structure of the QBX-discretized layer potential operators. In 
addition to being of mathematical interest, the spectral structure also plays an important role in determining 



the performance of iterative methods such as GMRES [Saad and Schultz 1986 when used to solve integral 
equations. 

We concentrate here on the double layer potential Dfi(x) which plays an important role, for example, in 
solving the Dirichlet problem for the Helmholtz equation (at a non-resonant frequency k) in the interior fl~ 
of r. Given Dirichlet data f{x), representing as a double layer potential 

(j){xo) = D^{xo) 



for xq ^ VL , and using the jump relation (18), we obtain the equation 



{-\+D)t,{x)^f{x) (24) 

for a; g r. 

On smooth boundaries, D is continuous and hence compact, with a discrete, bounded spectrum that has 



a unique accumulation point at zero Colton and Kress 1998 . In other words, it is a smoothing operator 



that damps out the high frequency modes in the density ^. We will show below that QBX is able to preserve 
all of these properties, most critically the spectral clustering at zero. 

We assume we have a grid on F with M panels and q points per panel, and that we consider a density 
that lives in the space of piecewise {q — l)th order polynomials over the M panels. We assume the double 
layer potential is computed using QBX, with values output on the same grid, corresponding to a discrete 
Mq X Mq matrix. 

Now, because the truncated p-term expansion represents locally smooth functions obeying a band limit 
related to p, high-frequency components of the density are either attenuated or aliased to lower frequencies 
as they transition from the source grid to the expansion. Empirically, attenuation is the dominant effect. 

This behavior has several consequences. A beneficial feature is that QBX responds very benignly to 
potentially erroneous high-frequency data that may be present in the discretized densities or geometries. 
Also, some spectral features are reproduced with no further effort. For example, when applied to the single 
layer potential, the QBX-based one-sided limit faithfully reproduces the spectrum of the continuous operator, 
accumulating at zero. 

Unfortunately, when computing the double-layer potential D using QBX based on the one-sided limit as 

-D/i,onc-sidcd K^) = lim Dhfi{x^) =F ^/^(a;) , 

high frequency components are attenuated in limj.±_j.2, _D;i/i(x*) but not, of course, in ^/i(a;). As a result, 
the spectrum of -Dh^onc-sidcd does not accumulate at zero. 



If one then solves the integral equation (24 1 iteratively, with I?/i,onc-sidcd computed in this manner, then 
an iterative method will converge rapidly up to the level of discretization error, at which point it will stall. 
Since one does not know a priori exactly what the discretization error will be, this is rather inconvenient. 



Fortunately, computing D using the two-sided averaging approach (23) as discussed in Section 3.2 yields 
a discrete operator with a spectrum accumulating at zero, because both limits are filtered. Matching this 
feature of the continuous operator allows iterative linear solvers converge to rapidly to solutions having 
residuals near machine precision even beyond the level of discretization error. Using two-sided averaging as 



in (23) may not be the only way to achieve this spectral behavior, but it is particularly convenient. 
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The preceding discussion applies only to the compact case. For hypersingular (finite-part) integrals or 
Hilbert-Riesz type operators such as i? in (22), there is no particular advantage in using the two-sided limit. 
Also, as discussed above, if jump conditions are not invoked, operators such as the (compact) single-layer 
potential can be represented faithfully by the one-sided procedure without difficulty. In the numerical results 
shown in Section [4j we have used two-sided averaging for all operators to which it applies. A more detailed 
discussion of the spectral properties of integral operators computed using QBX will be reported at a later 
date. 



4. Numerical experiments 

In this section, we illustrate the performance of QBX. We begin by describing some simple test geometries. 
We then present results for a variety of layer potentials, showing that high accuracy can be achieved even 
with modest-sized discretizations. Finally, we investigate the performance of QBX when used as part of an 
integral equation solver for a variety of Dirichlet and Neumann boundary value problems at various orders 
of accuracy. 

For the sake of convenience, we will denote the normal derivatives of the single and double layer potentials 

by 

S'a(x) := h{x) ■ \/ xSa{x), 
D'a{x) h{x) ■ VxDa{x). 

When x € r, the first is meant in the principal value sense and the second in the Hadamard finite-part sense. 



Four test geometries 

The four curves that we will use for our numerical tests are shown in Fig. |6] The ellipses of Figs. |6(a)[ 



6(b) 6(c) are given by 



-Yfn - ( cos(27ri) 

for Of = 1, 3, and 6, respectively, and the "starfish" of Figure [6 (d)| is given by 

/ sin(5 • 27ri)\ /cos(27rt)\ 



(25) 



(26) 



In each of these cases, t € [0, 1). 

We decompose the curves "f{t) = {x{t),y{t)) into panels, with x{t) and y{t) represented by a 16-term 
Legendre polynomial expansion. We generate an initial subdivision that is equispaced in t. To ensure the 
accuracy of the expansion, we sample the curve at 64 = 4 • 16 points per panel and compute Legendre 
expansion coefficients by Gauss-Legendre quadrature. 

Since we arc integrating with respect to the parameter t rather than arc length, we first determine 
whether the curve is well-resolved by studying the spectral decay of the Legendre coefficients of |7'(i)| 



V^ii^y + 72(^)^1 using the method of 



Klockner et al. 



2011 



We then determine the energy contained 



in the tail of the series (i.e. in modes 16 and above in our case). If the estimated residual exceeds 10"^^, the 
panel is bisected. A panel is also bisected if its length (as computed by integrating |7'|) is more than twice 
that of its neighboring panels, to avoid the issues described in Section [3.4. 1| 

A source oversampling factor (see Section 3.4) of 6 is used throughout, that is g — 6 • 16. A factor of 



2 is included to allow adjacent panels to differ in length by factors of two, and an additional factor of 3 
is included to ensure that the second term in the estimate ( 10 ) is negligible. More judicious oversampling 
strategies will be considered at a later date. 



4-. 2. Layer Potential Evaluation 

Our first set of tests examines the ability of QBX to compute a range of standard and non-standard layer 
potential operators to high precision. We consider the operators S, dxS, dyS, dxxS, dxyS, dyyS (target 
derivatives of the single-layer potential) as well as D, dxD, dyD (target derivatives of the double-layer 
potential). We also consider the layer potential induced by tangentially oriented dipoles (a source derivative 
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Figure 6. Test geometries and their panel subdivisions. 
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(a) Circle using K = 50 panels (b) 3-to-l ellipse using K = 50 panels 
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Table 1. Relative error of QBX integral operators in the and L°° norms compared to a high-accuracy 
reference solution, using local expansions of order p = 16. 
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(a) Sources and targets for the test of the solution 
of an exterior boundary value problem. 



(b) Sources and targets for the test of the solution 
of an interior boundary value problem. 



Figure T. Setup of the integral equation test cases, shown with the 'starfish' geometry of Figure 6(d) 
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Table 2. Convergence in the norm of the solution evaluated at a set of targets after solving a boundary value 
problem using an integral equation and QBX on the circle of Figure 6(a) GMRES iteration counts are shown in 
parentheses next to the error data. "EOC" is the empirical order of convergence, obtained by a log-least-squares 
fit of the shown errors. 



in the tangential direction), which we denoted earlier by by R. R is the analog for Helmholtz potentials of 
the Hilbert transform in two dimensions or the Riesz transform in three dimensions. 

We apply each of these operators to the density a{t) — sin(107r<) and compare the computed result 
to a reference solution in the and L°° norms. The Helmholtz parameter was chosen as fc = 0.5. The 
computations were carried out with local expansion order p = 16. We obtained our reference solution 
by using adaptive Gaussian quadrature with tolerance 10~^^ in quadruple precision with target points at 
distances 10^^, 10~^/2, and 10~^/4 from the curve along the normal on either side. We then computed 
one-sided limits w+ and v~ on each side by third-order Richardson extrapolation. We computed the value 
{v^ +v~)/2 as the reference solution for principal value or finite-part on-surface integrals. Results are shown 
in Table [T] confirming that high accuracy is achievable with modest computational effort, as expected from a 
rapidly convergent scheme. We further note that operators involving derivatives with tangential components 
to the curve are either hypersingular or bounded (but not compact). Since differentiation is ill-conditioned, 
one should expect some loss of accuracy with successively higher derivatives. 

4-.3. Integral equation solvers 

Our second and perhaps more important test examines the suitability of QBX in the context of solving 
integral equations of the second kind. For each interior domain, we define an exact solution as the field 
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Table 3. Convergence in the norm of the solution evaluated at a set of targets after solving a boundary 
value problem using an integral equation and QBX on the 3-to-l ellipse of Figure [6(b)[ GMRES iteration counts 
are shown in parentheses next to the error data. "EOC" is the empirical order of convergence, obtained by a 
log-least-squares fit of the shown errors. 
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Table 4. Convergence in the l'^ norm of the solution evaluated at a set of targets after solving a boundary 
value problem using an integral equation and QBX on the 6-to-l ellipse of Figure 6(c) GMRES iteration counts 
are shown in parentheses next to the error data. "EOC" is the empirical order of convergence, obtained by a 
log-least-squares fit of the shown errors. 
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Table 5. Convergence in the I norm of the solution evaluated at a set of targets after solving a boundary value 
problem using an integral equation and QBX on the "starfish" geometry of Figure 6(d) GMRES iteration 
counts are shown in parentheses next to the error data. "EOC" is the empirical order of convergence, obtained 
by a log-least-squares fit of the shown errors. 
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(b) Convergence in the P norm of the solution evaluated at a set of 
targets after solving a boundary value problem using an integral equation 
and QBX on the 'teardrop' geometry of Figure [8(a)| 



Figure 8. Integral equation tests on a 'teardrop' geometry with a corner. 
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induced by a collection of singular point sources in the exterior. For each exterior domain, an exact solution 
is constructed using singular point sources in the interior. Given the exact solution, we compute either 
Dirichlet or Neumann data and solve the corresponding boundary value problem using an integral equation. 
We then test the accuracy of the solution at a set of target points. Figure [7] illustrates the geometry of our 
tests. 

For the Dirichlet problem, we use the combined-field representation u = —Dk<J + aSkO' Colton and Kress 



1998 , which leads to the second-kind equation 



=F^cr + aSkCr - D^a = f, 



(27) 



for a, where / is the Dirichlet data obtained from the manufactured solution. Since the wave number is not 
large, we choose a = i throughout this section. The sign corresponds to the exterior problem and the 
(— ) sign to the interior problem. The subscript k in Sk and is used to emphasize that the underlying 
Green's function is that for the Helmholtz equation with Helmholtz parameter k. 

For the Neumann problem, we use a slight variation on the well-known combined-field representation 



Brakhage and Werner 1965 Bruno et al. 2012 Colton and Kress 1998 Leis 1965 Panic 1965 



u = SkO- - aDkSoa, 

where once again subscripts of k indicate the use of the Helmholtz kernel with parameter fc, and a subscript 
of indicates the use of a Laplace kernel. This representation leads to the second-kind integral equation 



1 



S'k<J - aD'kSoa = / 



(28) 



for (T, where / is the Neumann data obtained from the manufactured solution. Since QBX can integrate 



hypersingular kernels, we use (28) directly. One may also use the Calderon projection identity 

D'^So = -I /A -f S'^S'^ 



(29) 



Nedelec 2001 and some algebra to avoid hypersingular operators. 

Given our boundary discretization, we assume the unknowns are point values of a at the source nodes 
and enforce the integral equation at the same nodes, corresponding to a Nystrom method. We use GMRES 
to solve (27) or (28) iteratively and QBX to carry out the matrix- vector products. We set the GMRES 
tolerance to 10^^'Mndependent of the order of accuracy of the QBX-based quadrature. 



Following the work of Bremer 2012 , we use as unknowns the density values multiplied by the square root 
of the corresponding quadrature weight. This has the effect that the discrete P inner product approximates 
the continuous inner product and results in much improved conditioning, especially in the presence of 
widely varying panel sizes. (This is critical in geometries with corners, as discussed in the next section.) 

After solving for ct in (27) or (28), we use the corresponding representation to evaluate the potential u 



at a number of target points in f2. We then compare those values to the exact (manufactured) solution and 
compute the relative error. 

Results for the geometries described in Section 4T are shown in Tables [2j [3j [4j and [5] We observe that, 
while slightly more erratic, the results for the Neumann operator exhibit the loss of one order of accuracy, 
as expected since we have used QBX for a hypersingular kernel. 

Of particular note is the fact that, as predicted in Section [375} GMRES iteration reaches a residual of 
10~^* with a modest number of iterations even for low order accurate discretization, in nearly all cases. This 
makes QBX-based solvers particularly robust. 



4-. 3.1. Non-smooth geometries 

In the derivation of QBX, we have assumed that the layer potential is locally smooth, so that an expansion 
in Bessel functions is rapidly convergent. Since many engineering problems involve geometries with corners 
(and therefore potentially non-smooth densities and unbounded layer potentials), it is of interest to study 
whether QBX can be used effectively for such problems as well. 



Without knowing the precise singularity in the density, it is shown in [Bremer[ 2012 Helsing and Ojala 



2008bj that high-order polynomial approximation combined with high-order quadrature on a dyadically 



19 



refined mesh yields high-order accuracy. Thus, the only question is whether QBX can evaluate layer potentials 
on such structures without excessive work. To this end, we consider a 'teardrop' shape with a single corner, 
described by the parametrization 



7(t) = 1.7 



sin(7ri) — 0.5 
I cos(7rt)(7rt — 7r)7rt 



The curve and its discretization are shown in Figure 8(a) We have dyadically refined the boundary toward 
the corner until the smallest panel lengths are less than 10~^ on each side. Carrying out the same type 
of experiment as in the preceding section, we obtain the results in Table |8(b)[ The apparent drop in 
convergence order for p = 5 can be attributed to the error made in halting dyadic subdivision at e = 10^^. 
These experiments demonstrate that there are no significant obstacles to using QBX in this context. The 
method converges here precisely because the layer potential is locally smooth on finer and finer length scales 
as one approaches the corner. 



5. Generalizations and implementation issues 

Our goal in this paper has been to present a new approach to quadrature, which we refer to as QBX 
('quadrature by expansion'). While we have largely limited our attention to the Helmholtz equation in two 
dimensions, it should be clear that the overall approach is independent of dimension as well as the precise 
nature of the governing Green's function. In fact, the QBX approach is far more general, extending to 
kernels that are not directly connected to a partial differential equation. These extensions are discussed 
The method gives rise to a number of important and interesting questions regarding 



Klockner 2012 



efficiency, robustness, and automatic adaptivity. 

As for implementation, the main issue we have ignored here is computational cost. 
Section 



As presented in 



3.3 the asymptotic complexity of QBX is 0{NNt), where N is the number of source points, Nf is 



the number of target points. Neglecting numerous opportunities for optimization, a straight implementation 
of the algorithm of Section |3.3| can apply a single-layer operator at order p = 5 to a density on 7680 source 
nodes with 1280 targets in 0.6 seconds using 16 cores of a 2.93 GHz Intel Nehalem machine. Fortunately, 
QBX can be accelerated using the fast multipole method (FMM) or any other hierarchical fast algorithm 
Cheng et aL| |2006[ [Greengard and Rokhlinl |1987| . The cost is then 0{N\ogN + Nt log Nt ). This coupling. 



further cost savings, as well as extensions to three dimensions, are discussed in Greengard et al. 2012 



To give an indication of the achievable speedups, preliminary implementations show that the cost of an 
FMM-based QBX scheme for a layer potential is about twice that for a point-to-point FMM procedure. 
In particular, layer potentials with tens of thousands of discretization points are computed in seconds on 
a single CPU core. A variety of other optimizations are also possible: using direct evaluation for distant 
panel interactions and QBX for near neighbors only, adaptive oversampling to ensure accuracy of the local 
expansion coefficients with highly adaptive and irregular panel sizes, using the sample local expansion for 
several nearby target points, etc. 



6. Conclusions 

QBX permits the rapid, high-order accurate evaluation of layer potentials in a manner that is remarkably 
easy to implement. It is based on the fact that the induced potential is smooth in the exterior or interior 
domain. The scheme is equipped with a complete convergence theory. With minor modifications, QBX 
can evaluate layer potentials at off-surface points arbitrarily close to the boundary. Since it is an extension 
of the scheme developed in Barnett and Nguyen 2012 for precisely that purpose, this is not a surprise. 



QBX presents an opportunity to develop a robust set of software tools for evaluating integral operators with 
singular or weakly singular kernels, with application to a broad range of large-scale simulations in physics 
and engineering. 
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